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Abstract 

We present a constrained formulation of Dedner et al's hyperbolic/parabolic divergence cleaning scheme for 
enforcing the V • B = constraint in Smoothed Particle Magnetohydrodynamics (SPMHD) simulations. 
The constraint we impose is that energy removed must either be conserved or dissipated, such that the 
scheme is guaranteed to decrease the overall magnetic energy. This is shown to require use of conjugate 
numerical operators for evaluating V • B and Vijj in the SPMHD cleaning equations. The resulting scheme is 
shown to be stable at density jumps and free boundaries, in contrast to an earlier implementation by Price & 
Monaghan (2005). Optimal values of the damping parameter are found to be a = 0.2-0.3 in 2D and cr = 0.8- 
1.2 in 3D. With these parameters, our constrained Hamiltonian formulation is found to provide an effective 
means of enforcing the divergence constraint in SPMHD, typically maintaining average values of /i|V-B|/|B| 
to 0.1-1%, up to an order of magnitude better than artificial resistivity without the associated dissipation 
in the physical field. Furthermore, when applied to realistic, 3D simulations we find an improvement of 
up to two orders of magnitude in momentum conservation with a corresponding improvement in numerical 
stability at essentially zero additional computational expense. 



1. Introduction 

A key problem in numerical magnetohydrodynamics (MHD) is maintenance of the divergence constraint, 
V • B = from Maxwell's equations. If this is not maintained, a spurious force parallel to the magnetic 
field appears which can lead to numerical instability. A variety of methods have been developed to combat 
numerical divergence error, including Brackbill and Barnes [5] projection method, Evans and Hawley's [11] 
constrained transport, and Powell's [24l [25j eight wave approach or variants thereof. In general, these 
methods either aim to "clean" any divergence of the magnetic field that has been generated, or to alter 
the MHD formulation so that the divergence constraint is satisfied by construction. Toth [40] provides an 
excellent comparison of these schemes. 

However, even methods such as constrained transport which guarantee divergence free magnetic fields 
only do so in a particular discretisation. This means that numerical artefacts may still be present in different 
discretisations — such as those used in the force terms. The goal of all methods aimed at maintaining the 
divergence-free constraint is not to keep V • B exactly zero, but rather to prevent the growth of these 
numerical artefacts. 

Several attempts have been made to enforce the divergence constraint in the Smoothed Particle Hy- 
drodynamics (SPH) implementation of magnetohydrodynamics (SPMHD). One method is to formulate the 
magnetic field in terms of the Euler potentials [38] a and /3, setting B = Va x V/3, which guarantees zero 
divergence of the magnetic field by construction. Due to the Lagrangian nature of SPH, the scalar variables 
are advected exactly, which means the magnetic field can be reconstructed simply from the particle posi- 
tions relative to the initial conditions. While this method has found reasonable success (e.g. [13 [30]), the 
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range of problems it is suitable for is limited. For example, winding motions cannot be modelled past one 
rotation as the field is essentially "reset" with each turn. Price [28 investigated use of the vector potential 
formulation of the magnetic field as a way to overcome these limitations while still retaining the guarantee 
of zero physical divergence in the field. However, formulating the equations of motion in terms of the vector 
potential were found to be unstable, and there were significant difficulties were found with the time evolution 
of the vector potential, from which Price concluded that this was not a viable approach. 

Divergence cleaning schemes have also been used in SPMHD. The simplest approach is to rely on artificial 
resistivity to restrict growth of the divergence of the magnetic field. Artificial resistivity was introduced by 
Price and Monaghan [34j henceforth PM05]) to capture shocks and discontinuities in the magnetic field (see 
[27] for a general discussion on discontinuities in SPH), and corresponds to adding physical diffusion terms 
of the form |32| [29] 



dB 

~dt 



r]\/'B = T^V (V • B) - T^V X (V X B), (1) 



resist 



where r] is the resistivity parameter (in SPMHD r] ~ o^B'^sig^ where h is the smoothing length, '^sig is the 
maximum signal velocity and is a dimensionless parameter of order unity). This means that a form of 
divergence cleaning may be obtained as a byproduct at no additional computational expense. Biirzle et al 
[6] found this to be the case in star formation simulations, where maintaining the divergence constraint has 
proved challenging [31 . The caveat to relying on artificial resistivity to control divergence error is that it 
diffuses both physical and unphysical components of the magnetic field (Eq. [T]). Resistivity should not be 
increased in strength just to address divergence errors, as doing so will also weaken the physical field. 

Dedner et al.'s [8^ hyperbolic divergence cleaning scheme has found popular use in both Eulerian (ie, 
[19], [41^) and Lagrangian codes ([13^, [23^). To facilitate cleaning of divergence errors, an additional field 
ip is coupled to the magnetic field. The Dedner et al. scheme was originally adapted to SPMHD by |PMQ5" 



but was not adopted for two main reasons: i) the reduction in divergence error was relatively small (a factor 
of ^ 2-3) and ii) on certain test cases it was found that it could lead to an increase in the divergence error. 
As such, its use was not recommended [c.f. [29]. 

Our aim in this paper is to provide a formulation of hyperbolic divergence cleaning for SPMHD that is 
guaranteed to be stable and ensures a negative definite contribution to the magnetic energy. This means 
that the divergence cleaning is guaranteed to decrease the errors associated with non-zero divergence of the 
magnetic field, in turn leading to a method that is suitable for general use in SPMHD simulations. 

We begin with a review of the MHD equations and the divergence constraint for continuum and SPMHD 
systems (^. In ^ we discuss hyperbolic cleaning as part of the ideal MHD equations, and in ^3.2, define 
an energy term associated with the ip field. Using this energy term, we derive a new form for the ?/;-evolution 
equation which conserves total energy (^3.2.1). In ^ the discretisation of hyperbolic cleaning into SPMHD 
is discussed and we show how the constraint of energy conservation can be used to construct a formulation 
that is numerically stable. In particular, this leads to a requirement for the discretisation of V • B and S/t/j 
used in the induction and ?/^-evolution equations to form a conjugate pair (^4.2.1 and ^4.2.2). Importantly, 
we prove that the dissipative (parabolic) term in the evolution of t/j gives a negative definite contribution to 
magnetic energy (^4.3). Our new, constrained formulation of hyperbolic cleaning in SPMHD is then applied 
to a suite of test problems designed to evaluate all aspects of the algorithm and to derive parameter ranges 
suitable for general use (^. The final test (^5.7) is drawn from our current work on applying the method 
to star formation problems and shows that our technique performs well in practice, dramatically improving 
the accuracy and robustness of realistic SPMHD simulations in three dimensions. The results are discussed 
and summarised in ^ 



2. Smoothed particle magnetohydrodynamics 

2.1. Equations of ideal magnetohydrodynamics 

The ideal MHD equations are formulated by the coupling of the Euler equations of fluid flow with 
Maxwell's equations and the Lorentz force, under the assumption of the fluid being a perfect conductor with 
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zero net charge. Electric field contributions to the dynamics are therefore considered negligible, and the 
resulting equations of motion are given by the continuity, momentum, and induction equations. 



dt 
dv 

dt 
dB 

~dt 
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where p is the density, v is the fluid velocity, P is the thermal pressure, B is the magnetic field and d/dt is 
the comoving (Lagrangian) time derivative. The magnetic force is written as the divergence of the Maxwell 
stress tensor, 



V • M = —V • {-^BH + BB) = — [-V {^B^) + (B • V)B + B(V • B)] . 

flQ PO 



(5) 



The first term provides a magnetic pressure analogous to thermal pressure, while the second term gives the 
tension force perpendicular to the magnetic field lines. The third term occurs only for non-zero divergence 
of the magnetic field and results in an unphysical force parallel to the field lines [51 [TH [9] . 

2.2. The divergence constraint in MHD 

The divergence constraint enters the MHD equations only as an initial condition. 



d_ 
dt 



(V • B) = 0, 



or, with the "source term" approach [231 IE] corresponding to Eq. [ij 

d{V ■ B) 



dt 



■ V • (vV • B) = 0, 



(6) 



(7) 



which is identical in form to the continuity equation (Eq. [2| but with the density replaced by V • B. 

2.3. Discrete form of ideal MHD equations 

The discrete SPMHD representation of the ideal MHD equations are given bv [32ll34] 



Pa = ^rnbWab{ha), ha /^fa 
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(9) 
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where Wab{ha) = VF(|ro — r^l, /i^) is the SPH smoothing kernel, h is the smoothing length, m is the particle 
mass, the summation is over the particle's neighbours and we use the subscripts a and b to refer to the 
particle index. In the variable smoothing length formulation of SPH, ^1 is a term that arises due to gradients 
in the smoothing length and Eq. [S] is a non-linear set of equations that is solved iteratively for both h and 
p (for details see [33l l35]). 
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2.4- V • B in SPMHD 

There are two basic operators for calculating first derivatives in SPH, which we refer to as 'difference' 
and 'symmetric' [29 . For the divergence of the magnetic field, the difference operator is given by 



(V • B)^ = --^ Yl (Ba - Bb) • VaWabiha), 



(11) 



while the symmetric operator is 



(V-B)„ = Pa^mf, 

b 



^bPb 



■ "^aWabih) 



(12) 



2.5. Tensile instability correction 

Price and Monaghan [33] showed that the discrete SPMHD momentum equation (Eq. |9| can be derived 
self-consistently from a Lagrangian using the discrete form of the continuity and induction equations (i.e., 
Eqs.jsjand 10, respectively) as constraints, ensuring that this is a consistent set of equations and that energy 
and momentum are conserved exactly. However, the conservative form of the momentum equation (Eq. [9| 
inherently contains the unphysical force parallel to the field lines (c.f. Eq. |5|. The effect is that a numerical 
instability — the "tensile instability" — occurs when ^B^ > P, causing particles to unphysically clump 
together. This can be effectively countered using B0rve et al's [3 approach to subtract the B(V • B) source 
term from the magnetic force. 



dVg 
dt 



VB 



^aPl 



• yaWab{ha) 



^bpi 



(13) 



which is added to Eq. |9] Since the instability manifests only when ^B^ > P, B0rve et al. [T introduce 
an adjustable parameter /3, showing that it is sufficient to use f3 = ^ to correct the instability in the 
magnetic pressure-dominated regime. Indeed, recently [1 have recommended using = ^ for general 



SPMHD calculations. However, we find in this paper (^5.4.5) that using /3 < 1 can produce numerical 
artefacts (c.f. Fig. 12). We therefore strongly recommend using P = 1 and adopt this throughout the paper 

Note that with P 



unless otherwise specified 
equivalent to Powell's eight wave approach [24j . 



1 the induction and momentum equations are formally 



2.6. Dissipative terms in SPMHD 

Dissipative terms are added to SPMHD in the form of artificial viscosity, resistivity, and thermal con- 
duction in order to treat discontinuities and shocks. The form of artificial viscosity and thermal conduction 
used in this work is that of Monaghan's [20], derived by analogy with Riemann solvers. The form of artificial 
resistivity used is that of Price and Monaghan [34] , given by 



dt 



Pa Yl rUb^^iBa - Bb)r ■ y^Wab- 



(14) 



diss 5 Pab 

Each particle has as set individually, and is allowed to vary according in the range ^ [0, 1] according to 



OiB^a , /iVxB^I |V-B«| 
max 



(15) 



where r = ha/cFBV^i^ with gb = 0.1. Thus resistivity is only strong at large gradients in the magnetic field 
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3. Hyperbolic divergence cleaning 



3.1. Hyperbolic divergence cleaning for the MHD equations 

Hyperbolic divergence cleaning involves the introduction of a new scalar field, that is coupled to the 
magnetic field by a term appearing in the induction equation, 



dB\ 



and the field xf) evolves according to 



dV^ 



-4V.B 



In the comoving frame of the fiuid, Eq. [THjand 17 combine to produce a damped wave equation 



a"(V-B) 



4V'(V-B) 



1 d{V ■ B) 
r dt 



0. 



(16) 



(17) 
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The equation above shows that this approach spreads divergence of the magnetic field like a wave away from 
a source, diluting the initial divergence over a larger area, enabling the parabolic (diffusion) term, — V^/r, to 
act more effectively in reducing it to zero. The wave speed, c/j,, is chosen to be the fastest speed permissible 
by the time step, typically equal to the speed of the fast MHD wave. A key consideration is setting the 
damping strength correctly to achieve critical damping of the wave, which maximises the benefit of wave 
propagation without damping being too weak. Dedner et al. suggested using 1/r = ChCr where Cr = 0.18, 
though this is problematic as Cr is not a dimensionless quantity. Instead, PM05 define 



h ' 



(19) 



where h is the smoothing length and a is a dimensionless quantity specifying the damping strength. PM05 
found that optimal cleaning was obtained for a G [0.4, 0.8] in their tests. A similar form was also adopted 
by Mignone and Tzeferacos [19 in their Eulerian code, who suggested values a G [0, 1]. 

3.2. Energy associated with the field 

For later purposes it will be useful to define an energy term associated with the i/j field, (here defined 
as the energy per unit mass). Specifically, the energy should be defined such that, in the absence of damping 
terms, any change in magnetic energy should be balanced by a corresponding change in e^. This is not 
merely a book-keeping exercise, as it will enable us to construct a formulation of hyperbolic divergence 
cleaning in SPMHD that is guaranteed to be stable. 



If we consider the closed system of equations formed by equations [16] and [T7| the total energy of the 
system can be specified according to 



E 



_2/ioP 



pdV. 



Conservation of energy in this subsystem implies 



d^ 
~&t 



B_ ^dB\ 



/ - ( 



de^ 
dt 



pAV = 0, 



(20) 



(21) 



where we have used the fact that d(/)dV)/dt 
to the time derivative of ?/^, giving 



0. We assume that the time derivative of can be related 



B /dB\ dV^ 

Mop W^/^ '^dt 



pdV = 0, 



(22) 



where x is an unspecified variable to be determined. Using Eqs. [T6]and[T7|in the absence of damping gives 



Pop 



Integrating the first term by parts, we obtain 



pdV = 0. 



LPoP 



(V • B)pdV / V^B • ds = 0. 

Po 



(23) 



(24) 



We take the surface integral in equation [24] to be zero. If the bounding surface is taken to be at infinity, 
then this assumption is reasonable as it should be expected that the amplitude of a divergence wave would 
be diluted to zero at such a limit. For closed systems, it is not clear how the surface term should be treated. 
However, similar surface terms appear in the standard SPH formulation and are treated by the addition 
of diffusion terms to capture discontinuities [27]. For this reason we investigated adding an artificial ip- 
diffusion term to account for ^/^-discontinuities, but found no particular advantage to using this in practice 



(see Appendix A ). 

From (24) we conclude that energy conservation requires x 



energy of the field should be defined according to 



2popc% 



2 • 



tjj/jiiopcj^^ and therefore that the specific 



(25) 



3.2.1. Energy conservation as part of the ideal MHD equations 

Considering total energy conservation with hyperbolic divergence cleaning included as part of the set 
of ideal MHD equations, additional terms relating to dp/dt appear in the preceding analysis (along with 
kinetic and other energy terms). Any terms not involving do not need to be considered as they conserve 
energy together [see 33 , so energy conservation reduces to the condition 



/ 



B 



dtp 



popcl dt 2poplcl dt 



ip'^ dp 

272" 



pdV = 0. 



(26) 



The first two terms balance each other, however, the third term remains. There are several possible ap- 
proaches to ensuring total energy conservation with respect to this term. One approach which we explored 
was to derive the MHD+cleaning equations from a Lagrangian that includes the term. The result is 
that an additional isotropic pressure term, —^ip'^ /{pocjj, appears in the momentum equation. Since it is 
undesirable to change the physical forces in the system, we instead adopt a simpler approach, which is to 
slightly modify the evolution equation for ip. 

From the continuity equation (Eq. [2| , we can deduce that the third term in Eq. 26 will be balanced by 
replacing Eq. 17 with 

|^ = -4V.B-^-i^V.v. (27) 



4. Hyperbolic divergence cleaning in SPMHD 

4.1. Hyperbolic divergence cleaning in SPMHD 



Hyperbolic divergence cleaning in SPMHD can be constructed for either the difference (Eq. 11) or 
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While both measure 



symmetric (Eq. 12) measure of V • B by using the appropriate operator in Eq. 
the divergence of the magnetic field, they do not provide the same measurement. For example, if a random 
distribution of particles is given a uniform magnetic field, the difference form will measure precisely zero 
— since the magnetic field is equal for all particles — but the symmetric form will not because it will 
reflect the disordered particle arrangement. Thus, it may be expected that the difference operator in general 
gives a more accurate measure of V • B and should be the operator used for cleaning. On the other hand. 
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it is the symmetric form which is used in the momentum equation (Eq. [9| and correspondingly in the 
tensile instability correction (Eq. [T3|, and cleaning in this operator may be more effective at improving the 
conservation of momentum. Thus in the context of divergence cleaning, it is not clear a priori which of the 
two should be preferred. This is one of the questions we will seek to answer in our tests. 

It is also not clear how the operator for VV^ should be chosen. In |PMQ5 a difference operator was used 
for both V • B and Vijj. However, the choice of operator for VV^ turns out to be an important issue in 
ensuring a stable method. 

4-2. Energy conservation of discretised hyperbolic divergence cleaning 

The key constraint we wish to impose on our divergence cleaning scheme is that the total magnetic 
energy should never increase due to cleaning. That is, any magnetic energy transferred into the ip-fi.e\d 
should either be conserved or dissipated. Specifically, in the absence of damping terms, the propagation of 
divergence waves should conserve energy, not only in the continuum limit but also in the discrete system. 



We can thus use the derived in ^3.2 to derive stable formulations of hyperbolic divergence cleaning for 
SPMHD — for either difference or symmetric V • B operators. 



As in §3.2[ we first consider only the subsystem formed by Eqs. [16] and 17 This means that for the 



moment we do not consider additional terms related to dp/dt (these are discussed in ^4.2.3). The total 



energy of the subsystem (Eq. [20]) can be discretised by writing the integral as a sum and replacing the mass 



element pdV with the particle mass m, giving 



PoPa 



POpaCl 



Assuming that the total energy of the subsystem is conserved, we have 

B^, / dBa \ , V^a dllJa 



dE V- 

-dt^^"^^ 



PoPa 



dt 



POpaCl dt 



= 0. 



(28) 



(29) 



4-2.1. Hyperbolic cleaning with difference operator for V • B 

If we choose to clean using the difference operator for V • B, then the SPMHD version of Eq. 17 in the 
absence of the damping term is given by 



dt 



V rUb {Ba - B5) • VaWabiha). 
ilapr ^ 



Using this in Eq. (29), we have 

rrig ^ ^ ( dBg 
PoPa \ dt 



E 



Expanding the right hand side into two separate terms gives 

m„, „ f dB„, \ niamb 



= _ V — V'aTT^ V nib {Ba - Bfc) • VaWabiha). 



(30) 



(31) 



E 



B„ 



dt 



EE 



EEilf^^<^B,.V.M/.,(M, 



(32) 



where by swapping the arbitrary summation indices a and b in the second term on the right hand side and 
using the anti-symmetry of the kernel gradient (VaWab = —^bWba)^ we can simplify to find 



E 



POPa 



dt 



E 



POPa 




"^"^VaWabiha) + J^VaWab{hb) 



^bPb 



(33) 



This gives the SPMHD version of Eq. 16 in the form 



dB„ 

dt 



= -Pa XI 



^bPl 



^aWabih) 



(34) 



Thus, by choosing the difference operator for V • B, the symmetric operator for is imposed. That is, 
the total energy of the hyperbohc divergence cleaning scheme is only conserved if the operators are chosen 
to form a conjugate pair [c.f. [3 [28l |29]. This is an important improvement over the "PM05" implementation 
which used a difference operator for both. We demonstrate in ^that indeed the use of conjugate operators 
significantly improves the robustness and stability of our cleaning algorithm in practice. 

4.2.2. Hyperbolic cleaning with symmetric operator for V • B 

An energy-conserving formulation can also be constructed for divergence cleaning with the symmetric 



operator. That is, with Eq. 17 discretised according to 



dt 



b 



^aWabiha) 



B5 

^bpl 



^aWabih) 



the discrete version of Eq. [16] which must be used to conserve energy is constrained to be 

^dB " 



dt 



^aPa 



(35) 



(36) 



which again forms a conjugate pair. 



4.2.3. Hyperbolic cleaning as part of the SPMHD equations 



In ^3.2.1 



the evolution equation for was modified to include a — ^?/;(V • v) term (Eq. 27). This was 



done to conserve energy in the presence of dp/dt terms. The discretised form of V • v in Eq. 27 should 
therefore be the same as that used in the SPH continuity equation (see [2T), which leads to 



2ftaPa 



rribi^Ta - V5) • \/aWab{ha). 



(37) 



4.3. Energy loss due to damping 

For completeness, it is important to prove that the damping term in Eq. |27| will result in a negative 
definite energy change. Inserting the damping term into the total change of tjj energy, we see that 



dE\ 
~dtj 



damp 



E 



rrin 



POpaC% 



dt 



rria 



damp 



(38) 



which is indeed negative definite. These energy changes could be balanced with equivalent increases in 
thermal energy to keep the total energy constant. The issue with doing this is that the heat generated is not 
necessarily deposited in the same location as it was removed from the magnetic field, due to the transport 
of divergence errors inherent in the hyperbolic cleaning scheme. Thus, we do not add such heat gains as 
part of our method, although the term above can be used to keep track of the energy loss due to divergence 
cleaning. 



5. Tests 

We have designed our numerical tests to examine the following key aspects of our constrained hyperbolic 
divergence cleaning algorithm: 
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i) The importance of the energy-conserving, "constrained" formulation compared to a non-conservative, 
or "unconstrained", approach, 

ii) Whether or not cleaning using the symmetric V • B operator (^4.2.2) provides any advantage over use 
of the difference operator (^4.2.1), e.g. by improving momentum conservation. 



iii) Optimal parameter choices for a, 

iv) The practical effect of including the 



-^V^V • V term (Eq. 27). 



In particular, we have investigated these aspects both in isolation using simple idealised setups, as well as 
their combined effects in more realistic 2 and 3D simulations. Our goal is to verify the robustness of the 
algorithm for practical application in astrophysics, though it offers a general solution to maintaining the 
divergence constraint in SPMHD. 



As well as examining the divergence of the magnetic field using the operators given by Eq. 11 and 12 
we measure the divergence error in the standard manner for SPMHD with the dimensionless quantity. 



h\V-B\ 



(39) 



To prevent artificially high values where |B| ^0, a small parameter e is added to |B| in the denominator, 
where e ~ 1% of the maximum B-field value. We find this is only necessary for the Orszag-Tang vortex 
problem. 

All of the tests have been performed using a Leapfrog integrator with magnetic field integrated alongside 
the velocity and timesteps set according to the standard condition At < mina(Ccour^a/'^sig,o), where Ccour = 
0.2 and 

'^sig,a is the MHD fast wave speed on each particle. We therefore use cj^^ — maX(j(t'sig,a) the 



hyperbolic cleaning, except for the final test (^5.7) where Ch is individual to each particle. The damping 
parameter is chosen to be a = 0.4 (a = 0.8 for the final test), except for cases when it is varied to find 
optimal values. Unless otherwise indicated, we use the standard SPH cubic spline kernel for all tests with 
^fac = 1-2 in Eq. [8] corresponding to ~ 18 neighbours in 2D and 58 neighbours in 3D. The magnetic field 
is specified in units such that ji^ = 1 in the code [c.f. [32]. Artificial resistivity is only used where noted, in 



which case it is applied as described in ^2.6 



5.1. Divergence advection 

The simplest test we consider consists of divergence in the magnetic field artificially induced in the initial 
conditions and advected by a uniform flow. The test is performed in a two dimensional periodic domain 
with three dimensional magnetic and velocity fields (2.5D). The first version of this test is identical to the 



'divergence advection problem' from [8 , as generalised by PM05 We use this to illustrate the basic features 
of the hyperbolic/parabolic cleaning approach and to examine the optimal choice of a when the divergence 
error has a scale comparable to the numerical resolution. 

5.1.1. Setup 

The domain is a square area of fluid in the region x^y e [—0.5,1.5]. The system has uniform density 
p = 1, with pressure P = 6 and 7 = 5/3. The velocity field is v = [1, 1] and Bz = A perturbation 

is created in the x-component of the magnetic field of the form 



1 



(r/ro) -2(r/ro) 



1 



r < ro, 



(40) 



where r = \/ + and ro specifies the radial extent. We set up the problem using 50 x 50 particles on a 
square lattice, giving h = 1.2 Ax = 0.048. 
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Figure 1: A fluid with uniform velocity has a blob of divergence introduced to the initial conditions. In the top row, no cleaning 
is applied and the divergence blob is advected exactly with the flow. Undamped cleaning (purely hyperbolic) is applied to 
the centre row and the divergence in the magnetic fleld is spread through the system as a system of interacting waves. In 
the bottom row, damped cleaning (mixed hyperbolic/parabolic) is utilised and the divergence in the magnetic field is rapidly 
removed. 
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Figure 2: Average and maximum V • B in code units, measured using the difference operator (Eg. as a function of time 
for the divergence advection test with ro = Without cleaning, the divergence for this simple problem remains constant. 

Using undamped cleaning (purely hyperbolic), the maximum divergence is reduced with an increase in average throughout the 
system. With damped cleaning (mixed hyperbolic/parabolic), both average and maximum are rapidly reduced. 



5.1.2. Results 

Fig. [l] shows renderings of V • B at various times from three calculations: no cleaning, undamped 
cleaning (purely hyperbolic), and damped cleaning (mixed hyperbolic/parabolic) with ro = 1/a/8, following 
p . These three calculations illustrate the basic ideas behind the divergence cleaning scheme: In the absence 
of any cleaning (top), the magnetic field and its divergence perturbation is advected without change on the 
particles. With the addition of hyperbolic cleaning (middle), the divergence errors are spread in a wave-like 
manner throughout the domain. Finally, the addition of the parabolic damping term (bottom row) acts to 
rapidly diffuse the divergence error to zero. 

This is demonstrated more quantitatively in Fig. |2j which shows the average and maximum values of 
|V • B| as a function of time for the three calculations. While purely hyperbolic cleaning can be seen to 
quickly reduce the maximum divergence error, the average error increases. The parabolic damping means 
that both the average and maximum values are reduced by an order of magnitude in roughly one wave 
crossing time {t ~ 0.3), and by roughly 5 orders of magnitude after several crossing times (t > 2). After this 
time the divergence error continues to decrease, but at a much slower rate (this is more obvious in Fig. [s] for 
the case = h). We attribute the turnover in the decay rate to the rapid removal of the short wavelength 
errors by the cleaning scheme, leaving only slowly decaying long- wavelength modes. We have confirmed this 
interpretation by verifying that the transition to a slow decay is independent of timestepping, resolution 
and is similar using the quintic [29 instead of the cubic spline kernel. 

5.1.3. Optimal choice o f damping parameter in 2D 

As noted by the optimal choice of damping parameter, cr, for this problem with ro = l/>/8 is 

somewhat misleading, since in reality one expects divergence errors arising in simulations to have length 
scales of order the smoothing length. Thus, Fig. [3] shows the average and maximum V • B in a series of 
calculations employing ro = h and values of a between 0.1 and 0.6. The results are similar to those shown 
in Fig. |2j with best results obtained in this 2D case using a ~ 0.2-0.3. 

5.2. Static cleaning test: density jump 

Our second test is a variant on the divergence advection problem, with identical setup (ro = 1/ a/8) except 
that the right half of the domain has its density increased by a factor of two. The idea is to examine the 
refiection and refraction of the divergence waves as they transition between media of differing densities, as 
may frequently occur in applications of SPMHD. To simplify the test, we solve only the subset of equations 
given by Eqs. p!6}|T7| — that is, the system can only evolve due to divergence cleaning. 
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Figure 3: The effect of varying the damping parameter cr on the average and maximum V • B for the divergence advection test 
with ro = h. The best results for 2D are obtained for values between 0.2-0.3. 
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Figure 4: Results of the static cleaning test across a 2:1 density jump. Undamped non-conservative cleaning (top) increases the 
divergence of the magnetic field at the density jump, in turn leading to numerical instability (Fig. [5|. Using our constrained 
divergence cleaning method (bottom), the waves cross the density boundary without issue and the scheme remains stable. 
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Figure 5: Maximum values of V • B (difference) for the density jump test for the non-conservative formulation (left) and 
the new constrained divergence cleaning (right). The interaction between the divergence waves and the density jump for the 
non-conservative formulation is unstable, for both damped and undamped cleaning. Using constrained divergence cleaning is 
stable across the density jump, with damped cleaning reducing V • B as in previous tests. 



5.2.1. Setup 

The setup is performed in 2D with 25 x 50 particles on a square lattice in the left half of the domain 
{x < 0.5, p = 1), and 35 x 70 particles placed in the right half of the domain {x > 0.5, p = 2), with all 
particles of equal mass, giving a 2:1 density jump at x = 0.5. The actual density on the particles is found in 



the usual manner by iterating the smoothing length and density self-consistently as described in ^2.3 The 
velocity field is set to zero, all other system parameters are set as previously for the divergence advection 
test (^5.1[), and periodic boundary conditions are employed. 



5.2.2. Results 

Fig. [4] shows the propagation of purely hyperbolic (a = 0) divergence waves in this test using i) the 
non-energy conserving formulation with difference operators for both V • B (Eq. pT]) and V?/^ (Eq. [36|), and 
ii) our new constrained hyperbolic divergence cleaning scheme with a difference operator for V • B and the 
conjugate, symmetric operator for VV^ (Eq. [34]). The corresponding time evolution of the maximum |V • B| 
is shown in Fig. [5j Using the unconstrained formulation, the interaction of the divergence wave with the 
density jump causes amplification of the divergence errors (top row of Fig.|4|, in turn leading to exponential 
growth in the total energy and numerical instability (left panel of Fig.[5|. By contrast, our new conservative 
formulation remains stable and continues to reduce the divergence error throughout the domain (bottom 
row of Fig. [4] and right panel of Fig. |5|. 

5.3. Static cleaning test: free boundaries 

A further variant of the divergence advection test we consider replaces the periodic boundaries by a free 
boundary, since many applications of SPMHD involve free boundaries (e.g. the merger of two neutron stars 
[36], or studies of galaxy interactions [151 [16]). 

5.3.1. Setup 



The setup is identical to the divergence advection problem (^5.1 ) with ro = 1/a/8, except that the domain 
is a circular area of fluid with p = 1 for r < 1 and p = (no particles) for r > 1, set up using a total of 1976 
particles placed on a cubic lattice. The divergence perturbation is introduced at the centre of the circle, and 



the velocity field is set to zero. Rather than impose an external confining potential, we solve only Eqs. 16 -17 



without the full MHD equations, as in ^5.2 
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Figure 6: V • B of the static cleaning test using free boundaries. In the case of non-conservative cleaning (top row), the 
interaction of the divergence waves with the boundary cause unchecked divergence growth. Using constained cleaning (bottom 
row), the boundary interaction is not problematic. 



5.3.2. Results 

Fig. |6] shows the results of purely hyperbolic cleaning (a = 0) for this case. As in Fig. [4j the top 
row shows the unconstrained and non- conservative difference/difference formulation, while the bottom row 
shows results using the conservative difference/symmetric combination. Similar results are also found in 
this case, with divergence errors piling up at the free boundary in the non-conservative formulation leading 
to numerical instability, but our constrained formulation remaining stable. 

5.4' 2D Blast wave in a magnetised medium 

We now turn to tests that are more representative of the dynamics encountered in typical astrophysical 
simulations, beginning with a blast wave expanding in a magnetised medium. In this case the initial magnetic 
field is divergence-free, meaning that the only divergence errors are those created by numerical errors during 
the course of a simulation — rather than the artificial errors we have induced in the previous tests. Based 
on the results from the previous tests, in this and subsequent tests we apply cleaning only using constrained, 
energy-conserving formulations — that is, with conjugate operators for V • B and V?/^. We use this problem 
to the examine the effectiveness of the divergence cleaning in the presence of strong shocks, as well as to 
investigate whether cleaning should be performed using the difference or symmetric V • B operator. As with 
the divergence advection test, a key goal is to find optimal values for the damping parameter a. 

5.4-1- Setup 

The implementation of the blast wave follows that of Londrillo and Del Zanna [18]. The domain is a 
unit square with periodic boundaries, set up with 512 x 590 particles on a hexagonal lattice with p = 1. The 
fluid is at rest with magnetic field = 10. The pressure of the fiuid is set to P = 1, with 7 = 1.4, except 
a region of the centre of radius 0.125 has its pressure increased by a factor of 100 by increasing its thermal 
energy. An adiabatic equation of state is used. 
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Time Time 

Figure 8: Average and maximum of /i|V • B|/|B| as a function of time for the blast wave test. At t = 0.03, resistivity has 
reduced the average error by a factor of 4 compared to the control case, while divergence cleaning has reduced the average 
divergence error by a factor of 20. The maximum error has been reduced by a factor of 2 and 8, respectively. 
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Figure 9: V • B in the blast wave problem at t = 0.03 measured in code units using the symmetric V • B operator, showing the 
control case (left), V • B measured with the opposing operator to that used in the cleaning (centre) and V • B measured with 
the same operator used in the cleaning (right). Note in particular that the symmetric operator measures a divergence error 
around the leading edge of the fast MHD wave, even though the field is quite regular. 
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Figure 10: As in Fig. [o] but showing V • B measured using the difference operator. With this operator, no V • B is measured 
along the leading edge of the magnetic edge for the control and difference-cleaned cases. However, symmetric cleaning produces 
spurious divergence in this region when measured with the difference operator, because changes have been induced in the 
magnetic field to compensate for particle disorder. 

5.4.2. Results 

Fig. [t] shows the density and magnetic field lines at t = 0.03 for i) the control case without cleaning and 
no artificial resistivity (left), ii) including artificial resistivity (centre) and iii) no resistivity, but cleaned using 
the difference operator (right). At this time, the MHD fast shock has expanded to fill the domain, yet has 
not crossed the periodic boundaries to begin interacting with itself, and the three cases show only minimal 
differences in density structure. The average and maximum divergence error as a function of time are shown 
in Fig. [8] Although the density renderings at t = 0.03 are quite similar, we can see that adding divergence 
cleaning has reduced the average and maximum divergence error by a factor of 20 and 8, respectively at 
t = 0.03, compared to the control case, with factors of 5 and 4 improvement compared to the case with 
artificial resistivity alone. Thus, divergence cleaning is even more effective than resistivity at enforcing the 
divergence constraint. 

5.4-3. Operator choice for V • B 

To answer the question of whether there is any advantage to cleaning with the symmetric V • B operator, 
the blast wave problem was simulated for three cases: no cleaning; cleaning using the difference operator 
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Figure 11: Average and maximum h\V • B|/|B| for the blast wave test with varying damping strengths. The best results are 
obtained for values of a between 0.2-0.3. 



for V • B; and cleaning using the symmetric operator. The question is further comphcated by fact that the 
operator used for cleaning may differ from the operator used to measure the error. We therefore show V • B 
for these three cases measured with both the symmetric (Fig. [9| and difference (Fig. 10) operators, so that 
the effect of cleaning using one operator can be seen in both. 

The symmetric operator for V-B can be seen to pick up a non-zero divergence error on the leading edge of 
the magnetic wave from the blast (Fig. [9| despite the fact that the magnetic field shows no error in this region 



when measured with the difference operator (Fig. 10). This suggests that the symmetric operator is mainly 



reflecting the disordered particle arrangement. In turn, it can be seen that in this region, cleaning using the 
symmetric operator introduces divergence error when measured with the difference operator as it attempts 



to adjust the magnetic field based on the particle arrangement (centre panel of Fig. 10). Nevertheless, it 
is true that cleaning with the symmetric operator does produce the greatest reduction in the divergence 
when measured in the symmetric operator, and may still have potential advantages in terms of momentum 



conservation (this is examined further in ^5.5). However, we conclude that cleaning is best performed with 



the difference operator, since it shows not only the best results when measured with the difference operator 



(right panel of Fig. 10), but also an improvement even when measured with the symmetric operator (centre 
panel of Fig. [9|. 

5.4-4- Optimal damping values 

Fig. [TT] shows the average and maximum divergence error as a function of time for differing strengths 
of the damping parameter a in the range [0.1,0.6]. The best results are obtained with 0.2 < a < 0.3, in 
agreement with the other two dimensional tests. 

5.4-5. Tensile instability correction 

Finally, we noticed important consequences in this test concerning the /3B(V • B) correction for the 
tensile instability (§2^. Since using P = 0.5 is in principle sufficient to prevent the instability, its use has 
been suggested by [TJ[29]. However, we found this to be problematic in our simulations of the blast wave 
problem: Fig. 12 shows the density with overlaid magnetic field lines at t = 0.03 using P = 0.5, 0.75 and 1.0 
(left to right). With only /3 = 0.5 (left panel), irregularities can be seen to form in the densest parts of the 
Shockwave. These are not present when performing the full P = 1 subtraction (right panel). 

5.5. Orszag-Tang Vortex 

The final two dimensional test is the Orszag-Tang vortex [22j, which has been widely used as a test of 
MHD codes [e.g.[10][T2][39j. It consists of a magnetic vortex superimposed onto a velocity vortex generating 



17 



Figure 12: Density of the blast wave problem with overlaid magnetic field lines, without any divergence cleaning, but examining 
the impact of the —$'B(V • B) term used to correct the tensile instability. Results shown use $ = 0.5 (left), $ = 0.75 (centre) 
and (3 = 1.0 (right). Using only /5 = 0.5 is found to result in irregularities along the shock fronts, which are not present using 
= 1. Thus, using (3 = 0.5 is not recommended. 



several classes of interacting shock waves. The complex dynamics provides an excellent test of the constrained 
hyperbolic divergence cleaning method. To measure the effectiveness of the method in this case, the results 
are compared against that of simulations using artificial resistivity (with particle independent strengths 
as described in ^5.4.1) and Euler Potentials as measures of divergence control. This test is also used to 



examine whether or not cleaning using the symmetric operator for V • B provides any advantage in terms 
of momentum conservation. As previously, the damping parameter a is varied to find optimal values. 



5.5.1. Setup 

The problem is set up in a box with dimensions x^y G [0, 1] with periodic boundary conditions. The 
initial gas state is set to p = 25/(367r), P = 5/(127r), 7 = 5/3, with velocity field v = [— sin(27r?/), sin(27rx)]. 
The initial magnetic field is B = [— sin(27r?/), sin(47rx)]. All examples presented use 512 x 590 particles 
initially arranged on a hexagonal lattice. 



5.5.2. Results 



Fig. 13 shows the density (top), magnetic pressure (middle row), and V • B (bottom row) at t = 1.0 
for four cases: i) control, ii) using artificial resistivity, iii) employing Euler Potentials, and iv) applying 
divergence cleaning. This time is chosen because the divergence errors in the control case are large enough 
to produce small scale disturbances in the density and magnetic pressure fields. By adding resistivity or 
using Euler Potentials, the average /i|V-B|/|B| is decreased by an order of magnitude (c.f. second and third 



panels in bottom row of Fig. 13 and the left panel of Fig. 14). When divergence cleaning is used, the average 



divergence error is reduced by almost two orders of magnitude (red/dashed line in left panel of Fig. 14). In 



addition to the average and maximum divergence error for the above four cases. Fig. [M] also presents the 
results from a case where artificial resistivity has been applied in tandem with divergence cleaning. In this 
case, the average /i|V • B|/|B| is reduced by nearly an order of magnitude compared to resistivity alone, 
and when compared to the control case, this results in two orders of magnitude reduction in the average 
together with an order of magnitude reduction in the maximum. 



5. 5. 3. Cleaning using symmetric V • B 

Since the symmetric operator for V-B is used in the momentum equation and tensile instability correction 
term, it was hoped that its use for cleaning would confer some advantage over the difference measure by 



way of improved momentum conservation. However, as shown in Fig. 15, no significant difference in the 
momentum is found between cleaning with the symmetric operator compared to the difference operator. 



Fig. 16 shows the magnetic energy profile of the system for t < 0.5, where all test cases (control, resistivity. 



Euler Potentials, difference cleaning) yield the same profile, except for symmetric cleaning which shows a 
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Figure 13: The density (top row), magnetic pressure (middle row), and the difference measurement of V • B (bottom row) in 
the Orszag-Tang vortex at t = 1.0 comparing the control case (far left), including artificial resistivity (centre left), evolving the 
magnetic field using Euler Potentials (centre right), and applying the constrained divergence cleaning method (far right). 
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Figure 14: Average (left) and maximum (right) /i|V • B|/|B| in the Orszag-Tang vortex problem with (top to bottom in left 
panel) no divergence control; using Euler Potentials, adding an artificial resistivity, using divergence cleaning, and cleaning while 
including resistivity. Divergence cleaning has lower divergence error than when using Euler Potentials or artificial resistivity, 
and continues to reduce divergence error even when used in combination with artificial resistivity. 
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Figure 15: Total linear momentum for the Orszag-Tang 
vortex for divergence cleaning using the difference and 
symmetric operators of V • B. There is no significant 
distinction between the two. 



Figure 16: Magnetic energy as a function of time in the 
Orszag-Tang vortex test. Using the symmetric form of 
V • B for divergence cleaning leads to a 10% reduction 
in magnetic energy by t = 0.5 compared to the other 
schemes. 



^ 10% reduction in magnetic energy compared to the other solutions. This occurs due to the symmetric 
operator removing magnetic energy to compensate for irregularities in pa rticle position (which begin to 
occur at t ~ 0.15). Furthermore, although we have already shown in ^5.4.3 that use of /3 = ^ in the tensile 
instability correction could result in numerical artefacts in the blast wave test, we also found large errors in 
the density and magnetic field profiles when /3 = ^ is used in combination with symmetric cleaning on the 



Orszag-Tang problem. For these reasons, we recommend using (3 
difference V • B operator. 



1 and applying cleaning only with the 



5.5.4- Optimal damping values 

As with the previous tests, the damping parameter a was varied to find the best results (Fig.[l7|), which, 
as previously, were obtained for 0.2 < a < 0.3 for this 2D test. 



5.5.5. Resolution study 

Finally, the Orszag-Tang vortex test was performed for a series of increasing resolution: 128 x 148, 
256 X 296, 512 x 590, and 1024 x 1182 particles. Divergence cleaning, without resistivity, was used for all 
cases. The densities of these runs at t = 1.0 are shown in Fig. 18, along with results obtained using the 
Athena code ^39^ with 1024^ grid cells. In the largest resolution case, high density islands begin to form in 
the solution. These features are also exhibited in the results from the Athena code, and can be seen at lower 
resolutions in SPMHD when the Euler Potentials are used (see Fig. 13 for an example). Fig. 19 shows the 
average and maximum divergence error (left and right panels, respectively). Though the maximum error 
remains similar for all cases, the average is seen to decrease with increasing resolution. 



5.6. Three dimensional divergence advection 

We now turn to 3D tests, beginning with a three dimensional generalisation of the divergence advection 
problem. In particular, we wish to determine the optimal values for a when the divergence waves propagate 
in three dimensions rather than two. 



5.6.1. Setup 

The principle of the test remains similar to 2D versions, except a cubic volume of fiuid is used in the 
region x^y^z e [—0.5, 1.5]. The initial velocity field is extended to v = [1, 1, 1] to add drift in the z-direction. 
The magnetic field remains as previously, = I/a/Stt, with a spherical perturbation introduced to the 
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Figure 17: Average (left) and maximum (right) divergence error in the Orszag-Tang vortex problem, varying the damping 
parameter a. The best results are obtained with values ~ 0.2 — 0.3. 




Figure 18: Density of the Orszag-Tang vortex at resolutions of 128 x 148, 256 x 296, 512 x 590, and 1024 x 1182 particles 
(left to right), with comparison to results obtained using the Athena code for 1024^ grid cells (far right). As the resolution is 
increased, high density islands begin to form which is also observed in results from the Athena code. 
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Figure 19: Average (left) and maximum (right) divergence error in the Orszag-Tang vortex at resolutions of 128 x 148, 256 x 296, 
512 X 590, and 1024 x 1182 particles. The maximum divergence error remains similar for the different resolutions, but the 
average divergence error decreases with increasing resolution. 
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Figure 20: Average and maximum divergence error in the 3D advection test for varying strengths of the damping parameter, 
(J. The best results are obtained for a ~ 0.8-1.2. 



x-component of the field as given by Eq. 



40j except now using r = y^x^ + + z^. Tiie radial extent tq = h 
is chosen to mimic a divergence error at the resolution scale. The density and pressure remain unchanged, 
with p = 1, P = 6, and 7 = 5/3. The problem is set up on a cubic lattice with 50^ particles. 



5.6.2. Optimal values of the damping parameter 

This test was performed for a G [0.2, 1.2] with results of the average and maximum divergence given by 



Fig. 20 The optimal cleaning is obtained for a ~ 0.8-1.2, which differs from the optimal values obtained 
for the 2D tests of a ~ 0.2-0.3. This is attributed to the hyperbolic wave spreading spreading over a volume 
instead of an area, thus being more effective than in our 2D tests, and therefore requiring a higher value of 
(J to achieve critical damping. 



5. 7. Gravitational collapse of a magnetised molecular cloud core 

Our final test is drawn from our intended application: simulations of star formation that involve magnetic 
fields [37 . These simulations follow [30 , where an initial one solar mass sphere of gas with uniform magnetic 
field in the z-direction and in solid body rotation contracts under self-gravity to form a protostar with 
surrounding disc. However, at times near peak density, the magnetic field in the dense central region 
becomes strong and can produce high divergence errors. This has limited the range of initial magnetic field 
strengths which could be simulated, as if the divergence grows too large, the tensile instability correction 
term injects enough momentum into the system to erroneously eject the protostar out of its disc [31 . Thus, 
this simulation proves an excellent demonstration of the capabilities of the constrained hyperbolic divergence 
cleaning method to reduce divergence errors in realistic, 3D simulations. 

5.7.1. Setup 

The sphere of gas has radius = 4 x lO^^cm with uniform density p = 7.43 x 10~^^ g cm~^ and is set in 
solid body rotation with ft = 1.77 x 10~^^ rad s~^. A barotropic equation of state is used, as described in 
j3Q] . The magnetic field strength is set to give a mass-to-magnetic flux ratio of 5 times the critical value for 
magnetic fields to provide support against gravitational collapse. To avoid edge effects with the magnetic 
field, the sphere is embedded in a periodic box of length containing material surrounding the sphere set 
in pressure equilibrium with density ratio 1:30. This test uses only a minimal amount of resistivity, with 
aB ^ [0,0.1]. Self-gravity is simulated using a hierarchical binary tree where each node contains mutual 
nearest neighbours [2 , with gravitational force softening using the SPH kernel as described by [35 . The 
free fall time is ~ 24000 years. A sink particle is inserted once the gas density surpasses psink = 10~^^ g 
cm~^, and accretes particles within a radius of 6.7 AU. 
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Figure 21: Renderings of the column density of the star formation simulati on at t = 1.1 free fall times. The simulation without 
cleaning (left) suffers a dramatic loss of momentum conservation (c.f. Fig. |23| > induced by high divergence errors (c.f. Fig. |22[). 
By contrast, the simulation with our new divergence cleaning scheme applied (right) remains stable and launches a steady, 
collimated outflow [37) . 



5.7.2. Results 



Fig. 21 shows column density comparisons of simulations with (right) and without (left) divergence 
cleaning at t = 1.1 free fall time, showing that drastic improvements to the results are obtained by incor- 
porating divergence cleaning. Most importantly, the protostar remains stable in its disc. The average and 



maximum divergence error are both reduced by roughly an order of magnitude (Fig. ^2|, and this leads to 



a corresponding improvement in the momentum conservation of around two orders of magnitude (Fig. 23). 



5.7.3. Optimal sigma values 

This simulation was repeated for several values of the damping parameter in the range a G [0.2 — 1.2]. 
Optimal results were obtained for a ~ 0.8-1.2, which agrees with values found for the 3D advection test 
( ^5A2l ). 



5.7.4' Inclusion of the ^'^(V • v) term 

Adding -v) to the evolution equation for ip was motivated by energy conservation considerations, 

but the resulting question is what effect this has on divergence cleaning. The star formation simulation 
represents the ideal test case with which to examine this, with a large V • v present due to the gravitational 
collapse of the gas. We have performed this simulation both with and without this term, using a = 0.8, 
and found no distinguishable difference in the linear momentum, and average and maximum h\\/ • B|/|B| 
profiles. Similar results were obtained also found in the other tests. We conclude that, although this term is 
necessary for strict energy conservation, it has almost zero effect on the effectiveness of the cleaning scheme. 



6. Summary 

In this paper we have developed an implementation of Dedner et al's hyperbolic divergence cleaning for 
SPMHD that is constrained to be numerically stable and to always decrease the magnetic energy. To achieve 
this, we first defined the energy associated with the scalar tjj field (^3.2). This term was used to show that 
when the density varies over time, the evolution equation of ip should be modified to include a — • v) 

term in order to conserve energy. 

In [4.2 we derived an energy conserving formulation of divergence cleaning for SPMHD. By using the 
ip energy term, we showed that if a difference operator is chosen to discretise V • B in the dip/dt equation. 
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Figure 22: Average divergence error as a function of 
time for the star formation simulation, which shows that 
adding divergence cleaning reduces the divergence error 
by an order of magnitude. 



Figure 23: Magnitude of the total linear momentum in 
the star formation simulation. The system initially has 
zero net momentum, which increases due to the magnetic 
tensile instability correction and tree-based gravitational 
forces. After the protostar forms (t = 1), the momentum 
conservation in the divergence cleaning case is improved 
by two orders of magnitude over the control case. 
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Figure 24: Average and maximum divergence error in the star formation simulation, varying the damping parameter in the 
range a G [0.2, 1.2]. The best results are obtained with a ~ 0.8-1.2. 
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then the conjugate, symmetric operator for VV^ should be used ( p.2.1 ). Similarly, with symmetric V • B 



difference S/t/j should be used in the induction equati on (p.2.2 ). Use of conjugate operators was found to 



be the key to a numerical ly st able formulation. In §4.2.3[ we presented the correct SPMHD form of the 



-^?/^(V • v) term, and in ^4.3, demonstrated that parabolic damping will always lead to negative definite 
changes of energy. 

Tests of our constrained hyperbolic divergence cleaning were presented in ^ The selection of tests were 
for both 2 and 3D, and were designed to evaluate our method in isolation using simple, idealised systems 
and also in more realistic applications. Our idealised 2D tests consisted of a divergence advection test 



(^5.1[), and variants involving a density jump (^5.2) and free boundaries (^5.3). The more complex 2D tests 



were an MHD blast wave (^5.4) and the Orszag-Tang vortex (^5.5). A version of the divergence advection 



test extended to 3D was used in ^5.6 Results from the gravitational collapse of a molecular cloud core, 
representing our most challenging test case and a gauge of divergence cleaning applied to "real" applications, 
were presented in §5.7| From the results of these tests, we draw the following conclusions: 



Constrained hyperbolic/parabolic divergence cleaning provides an effective method of maintaining the 
divergence constraint in SPMHD, typically maintaining the average /i|V • B|/|B| to between 0.1-1%. 
The constrained formulation using conjugate operators for V • B and S/i/j is stable at density jumps 
and free boundaries, in contrast to previous implementations. 

We strongly recommend cleaning using the difference operator for V-B. Cleaning using the symmetric 
operator was not found to provide any advantage over the difference operator in terms of momentum 
conservation and was found to dissipate physical components of the magnetic field as well as the 
divergence error. 

Constrained divergence cleaning is more effective than artificial resistivity at reducing the divergence 
error, and still reduces the divergence error further when used in combination with resistivity. 
Divergence cleaning can provide an improvement of up to two orders of magnitude in momentum 
conservation when applied to realistic, 3D simulations. 

Optimal values for the damping parameter a were found to be cr = 0.2-0.3 in 2D and a = 0.8-1.2 in 
3D for all of the test problems considered in this paper. 

Addition of the —^i/jV-v term to the di/j/dt equation, while necessary for strict energy conservation of 
the hyperbolic cleaning equations, was found to have no noticeable effect, even in simulations where 
gas is strongly compacted. 

We found numerical artefacts in several problems when subtracting only — ^B(V-B) in the momentum 
equation to counteract the tensile instability. Instead, we strongly recommend using the full — B(V-B) 
correction. 

In summary, our constrained hyperbolic divergence cleaning scheme is a robust and effective method for 
enforcing the divergence constraint in SPMHD simulations, providing a pathway to accurate simulation of 
a wide range of magnetic phenomena in astrophysics and beyond. 
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Appendix A. Artificial ^-dissipation term 

Although the hyperbohc divergence cleaning method already includes a damping term to reduce V^, we 
have investigated the addition of a new dissipation term, analogous to artificial resistivity or viscosity, of 
the form 

=p,J2^,'-^i^P,-^,)F,,, (A.l) 

/ diss 5 Pab 

where \/Wab = ^abFab- This dissipation term is mainly designed to capture discontinuities in the tjj field. 



motivated by our neglect of the surface integral term in Eq. 24 The term is essentially an SPH expression 
for a diffusion term of the form r^^V^V^, where r]^ ex a^Chh, which in comparison to the damping term, acts 
more strongly to smooth relative differences in ?/^. This artificial ^/^-dissipation can be used in conjunction 
with the damping term, however since both the damping and diffusion terms dissipate ip^ it is important 
that values of and a be chosen carefully to avoid overdamping the system. For example, we found that 
in our two dimensional tests that propagation of divergence waves were damped too severely with = 1, 



and that using a^,cr = [0.1,0.2] or [0.2,0.1] yielded near critical damping (see Fig. A. 25). 
For this dissipation term, the energy loss is given by 



rrin 



/ J ^ 2 



This can be shown to be negative definite by splitting the RHS into two halves, performing a change of 
summation indices on the second half, then rejoining to obtain 



IE. 



oi^ sr^ {^a - ^bf ^ . 

I / =2 -^ab, (A.dj 

PoCh V Pab 



which, since Fab is negative for positive kernels, gives a negative definite contribution to the total energy 
(and conversely would give a positive definite heat contribution). 

Inclusion of the dissipation term was tried with all test cases presented in this paper. Similar reductions 
in the divergence error were obtained, however no results were improved beyond that of using the damping 



alone (Fig. A.25). 



27 



0.01 




Time Time 



Figure A. 25: Average and maximum divergence error when including the new, artificial dissipation term in the Orszag-Tang 
vortex test. Values of and a are chosen so that the combination is close to critical damping, however no benefit is noted 
over use of the regular damping term. 
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